% 老城街区搬迁补偿模型 - MATLAB实现
% 文件导入版本 - 修复版
% 日期：2025年

%% 主函数
function relocation_compensation_main()
    clc; clear; close all;
    
    % 数据导入
    fprintf('=== 老城街区搬迁补偿模型 ===\n');
    fprintf('正在导入数据文件...\n');
    
    % 设置文件路径
    plot_file = 'plot_data.xlsx';  % 地块信息文件
    coord_file = 'clicked_points.csv';  % 坐标文件
    
    try
        [plot_data, coordinates_data] = import_data_files(plot_file, coord_file);
        fprintf('数据导入成功！\n');
    catch ME
        fprintf('数据导入失败：%s\n', ME.message);
        fprintf('请检查文件路径和格式\n');
        return;
    end
    
    % 数据预处理
    [current_residents, available_plots] = preprocess_data(plot_data, coordinates_data);
    
    % 计算补偿方案
    fprintf('\n正在计算最优搬迁补偿方案...\n');
    [relocation_plans, total_cost] = calculate_relocation_plans(current_residents, available_plots);
    
    % 输出结果
    display_results(relocation_plans, total_cost, current_residents, available_plots);
    
    % 可视化结果
    visualize_results(relocation_plans, coordinates_data, plot_data);
    
    % 导出结果
    export_results(relocation_plans, current_residents, available_plots);
end

%% 数据导入函数
function [plot_data, coordinates_data] = import_data_files(plot_file, coord_file)
    
    % 导入地块信息文件
    fprintf('正在读取地块信息文件：%s\n', plot_file);
    
    [~, ~, ext] = fileparts(plot_file);
    
    if strcmpi(ext, '.xlsx') || strcmpi(ext, '.xls')
        plot_raw = readtable(plot_file, 'ReadVariableNames', true);
    elseif strcmpi(ext, '.csv')
        plot_raw = readtable(plot_file, 'ReadVariableNames', true);
    else
        error('不支持的文件格式：%s', ext);
    end
    
    % 提取数据并转换方位编码
    plot_data = zeros(height(plot_raw), 6);
    
    for i = 1:height(plot_raw)
        plot_data(i, 1) = plot_raw{i, 1};  % 地块ID
        plot_data(i, 2) = plot_raw{i, 2};  % 院落ID
        plot_data(i, 3) = plot_raw{i, 3};  % 地块面积
        plot_data(i, 4) = plot_raw{i, 4};  % 院落面积
        
        % 方位编码转换
        orientation_str = plot_raw{i, 5};
        if iscell(orientation_str)
            orientation_str = orientation_str{1};
        end
        plot_data(i, 5) = convert_orientation(orientation_str);
        
        plot_data(i, 6) = plot_raw{i, 6};  % 是否有住户
    end
    
    fprintf('地块信息读取完成：%d条记录\n', size(plot_data, 1));
    
    % 导入坐标文件
    fprintf('正在读取坐标文件：%s\n', coord_file);
    
    coord_raw = readtable(coord_file, 'ReadVariableNames', true);
    coordinates_data = table2array(coord_raw);
    
    fprintf('坐标信息读取完成：%d条记录\n', size(coordinates_data, 1));
    
    % 数据验证
    validate_data(plot_data, coordinates_data);
end

%% 方位编码转换函数
function code = convert_orientation(orientation_str)
    % 根据方位字符串转换为数字编码
    % 3=正南/正北（最优），2=东厢，1=西厢（最差）
    
    orientation_str = char(orientation_str);
    switch lower(strtrim(orientation_str))
        case {'南', 'south', '正南'}
            code = 3;
        case {'北', 'north', '正北'}
            code = 3;
        case {'东', 'east', '东厢'}
            code = 2;
        case {'西', 'west', '西厢'}
            code = 1;
        otherwise
            warning('未识别的方位：%s，默认设为东厢', orientation_str);
            code = 2;
    end
end

%% 数据验证函数
function validate_data(plot_data, coordinates_data)
    fprintf('\n数据验证中...\n');
    
    if size(plot_data, 2) ~= 6
        error('地块数据列数不正确，应为6列');
    end
    
    if size(coordinates_data, 2) ~= 3
        error('坐标数据列数不正确，应为3列');
    end
    
    residents_count = sum(plot_data(:, 6) == 1);
    available_count = sum(plot_data(:, 6) == 0);
    
    fprintf('数据验证通过：\n');
    fprintf('- 总地块数：%d\n', size(plot_data, 1));
    fprintf('- 有住户地块：%d\n', residents_count);
    fprintf('- 空置地块：%d\n', available_count);
    fprintf('- 坐标点数：%d\n', size(coordinates_data, 1));
end

%% 数据预处理函数
function [current_residents, available_plots] = preprocess_data(plot_data, coordinates_data)
    residents_mask = plot_data(:, 6) == 1;
    current_residents_data = plot_data(residents_mask, :);
    available_plots_data = plot_data(~residents_mask, :);
    
    % 处理当前居民数据
    current_residents = [];
    for i = 1:size(current_residents_data, 1)
        plot_id = current_residents_data(i, 1);
        
        coord_idx = find(coordinates_data(:, 1) == plot_id, 1);
        
        if ~isempty(coord_idx)
            x = coordinates_data(coord_idx, 2);
            y = coordinates_data(coord_idx, 3);
            
            is_streetside = is_plot_streetside(x, y, coordinates_data);
            is_convenient = is_plot_convenient(x, y, coordinates_data);
            
            current_residents = [current_residents; 
                plot_id, current_residents_data(i, 2), current_residents_data(i, 3), ...
                current_residents_data(i, 5), is_streetside, is_convenient, x, y];
        else
            current_residents = [current_residents; 
                plot_id, current_residents_data(i, 2), current_residents_data(i, 3), ...
                current_residents_data(i, 5), 0, 0, 0, 0];
        end
    end
    
    % 处理可用地块数据
    available_plots = [];
    for i = 1:size(available_plots_data, 1)
        plot_id = available_plots_data(i, 1);
        
        coord_idx = find(coordinates_data(:, 1) == plot_id, 1);
        
        if ~isempty(coord_idx)
            x = coordinates_data(coord_idx, 2);
            y = coordinates_data(coord_idx, 3);
            
            is_streetside = is_plot_streetside(x, y, coordinates_data);
            is_convenient = is_plot_convenient(x, y, coordinates_data);
            
            available_plots = [available_plots; 
                plot_id, available_plots_data(i, 2), available_plots_data(i, 3), ...
                available_plots_data(i, 5), is_streetside, is_convenient, x, y];
        else
            available_plots = [available_plots; 
                plot_id, available_plots_data(i, 2), available_plots_data(i, 3), ...
                available_plots_data(i, 5), 0, 0, 0, 0];
        end
    end
    
    fprintf('\n数据预处理完成：\n');
    fprintf('- 当前居民数量：%d户\n', size(current_residents, 1));
    fprintf('- 可用地块数量：%d个\n', size(available_plots, 1));
end

%% 判断是否为临街地块
function is_streetside = is_plot_streetside(x, y, coordinates_data)
    if x == 0 && y == 0
        is_streetside = false;
        return;
    end
    
    min_x = min(coordinates_data(:, 2));
    max_x = max(coordinates_data(:, 2));
    min_y = min(coordinates_data(:, 3));
    max_y = max(coordinates_data(:, 3));
    
    threshold = 100;
    is_streetside = (x - min_x < threshold) || (max_x - x < threshold) || ...
                   (y - min_y < threshold) || (max_y - y < threshold);
end

%% 判断交通便利性
function is_convenient = is_plot_convenient(x, y, coordinates_data)
    if x == 0 && y == 0
        is_convenient = false;
        return;
    end
    
    center_x = mean(coordinates_data(:, 2));
    center_y = mean(coordinates_data(:, 3));
    distance = sqrt((x - center_x)^2 + (y - center_y)^2);
    
    max_distance = max(sqrt((coordinates_data(:, 2) - center_x).^2 + ...
                           (coordinates_data(:, 3) - center_y).^2));
    
    is_convenient = distance < 0.4 * max_distance;
end

%% 计算搬迁方案主函数
function [relocation_plans, total_cost] = calculate_relocation_plans(current_residents, available_plots)
    n_residents = size(current_residents, 1);
    
    relocation_plans = [];
    total_cost = 0;
    
    for i = 1:n_residents
        resident = current_residents(i, :);
        fprintf('\n计算地块%d（面积%.0f平米，朝向编码%d）的搬迁方案...\n', ...
                resident(1), resident(3), resident(4));
        
        [feasible_plans, costs] = find_feasible_plans(resident, available_plots);
        
        if ~isempty(feasible_plans)
            best_plan_idx = select_best_plan(resident, feasible_plans, costs);
            best_plan = feasible_plans(best_plan_idx, :);
            best_cost = costs(best_plan_idx);
            
            relocation_plans = [relocation_plans; [resident(1), best_plan(1), best_cost]];
            total_cost = total_cost + best_cost;
            
            fprintf('  最优方案：搬到地块%d（面积%.0f平米），成本%.2f万元\n', ...
                    best_plan(1), best_plan(3), best_cost);
        else
            fprintf('  无可行搬迁方案\n');
        end
    end
end

%% 寻找可行搬迁方案
function [feasible_plans, costs] = find_feasible_plans(resident, available_plots)
    current_area = resident(3);
    current_orientation = resident(4);
    
    feasible_plans = [];
    costs = [];
    
    for i = 1:size(available_plots, 1)
        plot = available_plots(i, :);
        new_area = plot(3);
        new_orientation = plot(4);
        
        % 约束检查
        area_ok = (new_area >= current_area) && (new_area <= 1.3 * current_area);
        light_ok = new_orientation >= current_orientation;
        special_ok = ~((current_orientation == 3) && (new_orientation == 1));
        repair_area_ok = new_area >= 0.9 * current_area;
        
        if (area_ok && light_ok && special_ok) || ...
           (~light_ok && repair_area_ok && special_ok)
            feasible_plans = [feasible_plans; plot];
            costs = [costs; calculate_relocation_cost(resident, plot)];
        end
    end
end

%% 计算搬迁成本
function cost = calculate_relocation_cost(resident, new_plot)
    communication_cost = 5; % 万元
    
    current_area = resident(3);
    current_orientation = resident(4);
    new_area = new_plot(3);
    new_orientation = new_plot(4);
    
    % 面积损失成本
    area_loss_cost = 0;
    if new_area > current_area
        extra_area = new_area - current_area;
        daily_rent = get_daily_rent(new_orientation);
        area_loss_cost = extra_area * daily_rent * 365 / 10000;
    end
    
    % 修缮补偿成本
    repair_cost = 0;
    if new_area < current_area
        area_diff_percent = (current_area - new_area) / current_area * 100;
        repair_cost = repair_cost + area_diff_percent * 0.15;
    end
    if new_orientation < current_orientation
        repair_cost = repair_cost + 2.5;
    end
    repair_cost = min(repair_cost, 5);
    
    % 距离成本
    distance_cost = 0;
    if resident(7) > 0 && resident(8) > 0 && new_plot(7) > 0 && new_plot(8) > 0
        distance = sqrt((resident(7) - new_plot(7))^2 + (resident(8) - new_plot(8))^2);
        distance_cost = distance / 100 * 0.2;
    end
    
    cost = communication_cost + area_loss_cost + repair_cost + distance_cost;
end

%% 获取日租金
function daily_rent = get_daily_rent(orientation)
    switch orientation
        case 3
            daily_rent = 45;
        case 2
            daily_rent = 24;
        case 1
            daily_rent = 24;
        otherwise
            daily_rent = 24;
    end
end

%% 选择最优方案
function best_idx = select_best_plan(resident, feasible_plans, costs)
    if length(costs) == 1
        best_idx = 1;
        return;
    end
    
    [~, best_idx] = min(costs);
end

%% 显示结果
function display_results(relocation_plans, total_cost, current_residents, available_plots)
    fprintf('\n=== 搬迁补偿方案结果 ===\n');
    fprintf('总搬迁成本：%.2f万元\n', total_cost);
    fprintf('可搬迁居民：%d户 / %d户\n', size(relocation_plans, 1), size(current_residents, 1));
    
    if ~isempty(relocation_plans)
        fprintf('\n详细搬迁方案：\n');
        fprintf('原地块ID\t新地块ID\t成本(万元)\n');
        fprintf('--------\t--------\t----------\n');
        for i = 1:size(relocation_plans, 1)
            fprintf('%d\t\t%d\t\t%.2f\n', relocation_plans(i, 1), ...
                    relocation_plans(i, 2), relocation_plans(i, 3));
        end
    end
end

%% 可视化结果
function visualize_results(relocation_plans, coordinates_data, plot_data)
    if isempty(relocation_plans) || size(coordinates_data, 1) < 2
        fprintf('数据不足，跳过可视化\n');
        return;
    end
    
    figure('Name', '搬迁方案可视化', 'Position', [100, 100, 1200, 800]);
    
    % 绘制所有地块位置
    subplot(2, 2, 1);
    scatter(coordinates_data(:, 2), coordinates_data(:, 3), 50, 'b', 'filled');
    title('地块分布图');
    xlabel('X坐标');
    ylabel('Y坐标');
    grid on;
    
    % 区分有住户和空置地块
    subplot(2, 2, 2);
    residents_mask = plot_data(:, 6) == 1;
    resident_coords = [];
    vacant_coords = [];
    
    for i = 1:size(plot_data, 1)
        coord_idx = find(coordinates_data(:, 1) == plot_data(i, 1), 1);
        if ~isempty(coord_idx)
            if residents_mask(i)
                resident_coords = [resident_coords; coordinates_data(coord_idx, 2:3)];
            else
                vacant_coords = [vacant_coords; coordinates_data(coord_idx, 2:3)];
            end
        end
    end
    
    if ~isempty(resident_coords)
        scatter(resident_coords(:, 1), resident_coords(:, 2), 60, 'r', 'filled', 'DisplayName', '有住户');
    end
    hold on;
    if ~isempty(vacant_coords)
        scatter(vacant_coords(:, 1), vacant_coords(:, 2), 40, 'g', 'filled', 'DisplayName', '空置');
    end
    title('住户分布');
    xlabel('X坐标');
    ylabel('Y坐标');
    legend;
    grid on;
    
    % 搬迁成本分布
    subplot(2, 2, 3);
    if ~isempty(relocation_plans)
        bar(relocation_plans(:, 3));
        title('搬迁成本分布');
        xlabel('方案编号');
        ylabel('成本 (万元)');
        grid on;
    end
    
    % 搬迁流向图
    subplot(2, 2, 4);
    if ~isempty(relocation_plans) && size(relocation_plans, 1) > 0
        for i = 1:size(relocation_plans, 1)
            from_id = relocation_plans(i, 1);
            to_id = relocation_plans(i, 2);
            
            from_idx = find(coordinates_data(:, 1) == from_id, 1);
            to_idx = find(coordinates_data(:, 1) == to_id, 1);
            
            if ~isempty(from_idx) && ~isempty(to_idx)
                from_pos = coordinates_data(from_idx, 2:3);
                to_pos = coordinates_data(to_idx, 2:3);
                
                plot([from_pos(1), to_pos(1)], [from_pos(2), to_pos(2)], 'r-', 'LineWidth', 2);
                hold on;
                scatter(from_pos(1), from_pos(2), 60, 'r', 'filled');
                scatter(to_pos(1), to_pos(2), 60, 'g', 'filled');
            end
        end
        title('搬迁流向图');
        xlabel('X坐标');
        ylabel('Y坐标');
        grid on;
    end
    
    fprintf('\n可视化图表已生成\n');
end

%% 导出结果
function export_results(relocation_plans, current_residents, available_plots)
    if ~isempty(relocation_plans)
        result_table = table();
        result_table.OriginalPlotID = relocation_plans(:, 1);
        result_table.NewPlotID = relocation_plans(:, 2);
        result_table.RelocationCost = relocation_plans(:, 3);
        
        filename = sprintf('relocation_results_%s.xlsx', datestr(now, 'yyyymmdd_HHMMSS'));
        writetable(result_table, filename);
        fprintf('结果已导出到文件：%s\n', filename);
    end
end

%% 创建示例数据（用于测试）
function create_sample_data()
    % 创建示例地块数据
    plot_data = table();
    plot_data.PlotID = (1:20)';
    plot_data.YardID = (1:20)';
    plot_data.PlotArea = 80 + rand(20,1) * 40;  % 80-120平米
    plot_data.YardArea = plot_data.PlotArea * 0.8;  % 院落面积为地块面积的80%
    plot_data.Orientation = repmat({'南'; '东'; '西'; '北'}, 5, 1);
    plot_data.HasResident = [ones(10,1); zeros(10,1)];  % 前10个有住户，后10个空置
    
    writetable(plot_data, 'plot_data.xlsx');
    
    % 创建示例坐标数据
    coord_data = table();
    coord_data.PlotID = (1:20)';
    coord_data.X = rand(20,1) * 1000;  % 随机X坐标
    coord_data.Y = rand(20,1) * 1000;  % 随机Y坐标
    
    writetable(coord_data, 'clicked_points.csv');
    
    fprintf('示例数据文件已创建：plot_data.xlsx 和 clicked_points.csv\n');
end